
#ifndef incl_material_utilities_h
#define incl_material_utilities_h

#include <iostream>
#include <cmath>
#include <Eigen/Dense>

using Eigen::VectorXd;
using Eigen::MatrixXd;


inline  double  KronDelta(int a, int b)
{
  return ( (a==b) ? 1.0 : 0.0 );
}


inline  double  determinant3b3(VectorXd& F)
{
    return ( F(0)*(F(4)*F(8)-F(5)*F(7)) - F(3)*(F(1)*F(8)-F(7) *F(2)) + F(6)*(F(1)*F(5)-F(2)*F(4)) );
}

void   matrix2vector(const MatrixXd& mat,  VectorXd& vec);


void   vector2matrix(const VectorXd& vec, MatrixXd& mat);


void  adjust_deformationgradient_2D(int sss, bool finite, MatrixXd&  F);


void  adjust_deformationgradient_2D_Fbar(int sss, bool finite, MatrixXd&  F, double detFc);


void  get_smallstraintensor(const MatrixXd& F, MatrixXd& eps);


void  get_deviatoricpart(const MatrixXd& F, MatrixXd& eps);


void  defgrad2eps(VectorXd&  F,  VectorXd&  eps);


//  b = F*Ft
void  computeLeftCGtensor(VectorXd& F, VectorXd& b);


//  C = Ft*F
void  computeRightCGtensor(VectorXd& F, VectorXd& C);


void getAmtx2D(bool axsy, double cc[4][4], double stre[4], double aa[5][5]);


void getQmtx2D(bool axsy, double aa[5][5], double stre[4], double q[5]);


int  getMaterialID(std::string& matkey);



inline void Idev2D(double  Idev[][4])
{
   double  r1d3 = 1.0/3.0, r2d3 = 2.0*r1d3;

   Idev[0][0] =  r2d3;     Idev[0][1] = -r1d3;	 Idev[0][2] = -r1d3;    Idev[0][3] = 0.0;
   Idev[1][0] = -r1d3;	   Idev[1][1] =  r2d3;	 Idev[1][2] = -r1d3;    Idev[1][3] = 0.0;
   Idev[2][0] = -r1d3; 	   Idev[2][1] = -r1d3;	 Idev[2][2] =  r2d3;    Idev[2][3] = 0.0;
   Idev[3][0] =  0.0;      Idev[3][1] =  0.0;    Idev[3][2] =  0.0;	Idev[3][3] = 1.0;

 return;
}


inline void Idev3D(double  Idev[][6])
{
   double  r1d3 = 1.0/3.0, r2d3 = 2.0*r1d3;

   Idev[0][0] =  r2d3;     Idev[0][1] = -r1d3;	 Idev[0][2] = -r1d3;    Idev[0][3] = 0.0;    Idev[0][4] = 0.0;    Idev[0][5] = 0.0;
   Idev[1][0] = -r1d3;	   Idev[1][1] =  r2d3;	 Idev[1][2] = -r1d3;    Idev[1][3] = 0.0;    Idev[1][4] = 0.0;    Idev[1][5] = 0.0;
   Idev[2][0] = -r1d3; 	   Idev[2][1] = -r1d3;	 Idev[2][2] =  r2d3;    Idev[2][3] = 0.0;    Idev[2][4] = 0.0;    Idev[2][5] = 0.0;
   Idev[3][0] =  0.0;      Idev[3][1] =  0.0;    Idev[3][2] =  0.0;   	Idev[3][3] = 1.0;    Idev[3][4] = 0.0;    Idev[3][5] = 0.0;
   Idev[4][0] =  0.0;      Idev[4][1] =  0.0;    Idev[4][2] =  0.0;   	Idev[4][3] = 0.0;    Idev[4][4] = 1.0;    Idev[4][5] = 0.0;
   Idev[5][0] =  0.0;      Idev[5][1] =  0.0;    Idev[5][2] =  0.0;   	Idev[5][3] = 0.0;    Idev[5][4] = 0.0;    Idev[5][5] = 1.0;

 return;
}


#endif
